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Abstract. A first-principles description of nuclear systems along the drip lines 
presents a substantial theoretical and computational challenge. In this paper, we 
discuss the nuclear theory roadmap, some of the key theoretical approaches, and 
present selected results with a focus on long isotopic chains. An important conclusion, 
which consistently emerges from these theoretical analyses, is that three-nucleon forces 
are crucial for both global nuclear properties and detailed nuclear structure, and that 
many-body correlations due to the coupling to the particle continuum are essential as 
one approaches particle drip lines. In the quest for a comprehensive nuclear theory, 
high performance computing plays a key role. 
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1. Introduction 

To understand why nucleonic matter is stable is one of the overarching aims and 
intellectual challenges of basic research in nuclear physics. To relate the existence 
and properties of nuclei to the underlying fundamental forces and degrees of freedom, 
is central to present and planned rare isotope facilities, see for example Refs. [IHH]. 
Important properties of nuclear systems are binding energies, radii, density distributions 
of nucleons, spectra, and decays. These quantities convey important information on the 
individual-nucleon motion manifesting itself in the shell structure of nuclei, including 
the appearance and disappearance of magic numbers, interplay between high-j unique- 
parity orbits and natural-parity states, and the rapid changes of nuclear properties 
around the reaction thresholds. 

To relate the stability of nucleonic matter at various energies and length scaled 
to the underlying fundamental forces is a major quest for theoretical modeling. A 
multiscale approach is required which incorporates different degrees of freedom at 
relevant length and energy scales. Unfortunately, quantum chromodynamics (QCD), 
the underlying theory of the strong interaction, is highly non-perturbative in the energy 
region characteristic of nuclear structure. This region is governed by nucleonic (and 
sometimes mesonic) degrees of freedom. This requires an effective theory, usually framed 
as chiral effective field theory (EFT), which is consistent with the symmetries of low- 
energy QCD and the separation of scales relevant to the low-energy nuclear many-body 
problem. Linking different scales is far from easy. A key challenge is to understand 
the link between QCD and effective field theories. Interactions derived from these low- 
energy theories carry also a dependence on an energy scale, defined in terms of an energy 
cutoff A, that separates the Hilbert space of interest from its higher-energy complement. 
The cutoff is usually chosen so that it is possible to reproduce nucleon-nucleon (NN) 
scattering phase shifts up to ~ 300 MeV laboratory kinetic energy, which requires 
A ~ 500 — 700 MeV. For details, we refer the reader to the recent reviews [9l[T0] . 

An important question is to understand the role of many-nucleon forces from chiral 
EFT in the nuclear medium. Furthermore, using EFT-based interactions in a many- 
body environment entails the development of proper many-body theories that allow 
for first principle calculations. There are other issues as well, such as the estimation 
of theoretical errors. For example, most many-body methods applied to the nuclear 
many-body problem involve basis set expansions. The errors which arise due to basis 
truncations need proper quantifications and clarifications [TTHT^ . Finally, a proper 
link between first-principle methods (which are of limited use when very many degrees 
of freedom are at play) and approaches based on the density functional theory, is 
essential if one wishes to understand nuclei and nuclear matter from a bottom-up 
perspective p3|fT6]. 

f A simultaneous description of nuclei and neutrons stars in terms of hadronic degrees of freedom 
constitutes a daunting theoretical challenge as the length scale spans over 19 orders of magnitude, from 
several 10^^^ m (nuclear radii) to approximately 12 kilometers (neutron star radii). 
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Neutron-rich nuclei are particularly interesting for this endeavor. As the neutron 
number increases within a particular chain of isotopes, one eventually reaches the limit 
of stability, the neutron drip line, where isotopes with additional neutrons are not 
neutron-bound anymore. The appearance or absence of magic numbers, formation of 
neutron skins and halos, and detailed tests of shell structure at the limit of neutron-to- 
proton asymmetry can be probed via investigations of masses, radii, and excited states 
of neutron-rich nuclei. Examples of recent progress in the radioactive nuclear beam 
(RNB) science are measurements of masses with Penning traps coupled to radioactive 
beams and nuclear charge radii and moments using laser spectroscopy techniques. New 
experimental data in the neutron-rich territory are crucial for constraining theoretical 
models of nuclei and astrophysical processes, see, e.g., Refs. |T7l[IH]- To put things 
in perspective. Fig. [T] shows the expected experimental information on the Ca isotopic 
chain that will be obtained at the future Facility for Rare Isotope Beams (FRIB). 
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Figure 1. Expected experimental information on the calcium isotopes that can be 
obtained at FRIB. The limits for detailed spectroscopic information are around ^ 60. 
(Courtesy of Brad Sherill [19].) 

Having access to precise measurements of masses, radii, and electromagnetic 
moments for a wide range of nuclei allows to study trends with varying neutron excess. 
A quantitative description of various experimental data with quantified uncertainty 
still remains a major challenge for nuclear structure theory. Global theoretical studies 
of isotopic chains, such as the Ca chain in Fig. [H make it possible to test systematic 
properties of effective nuclear interactions. Such calculations also provide critical tests of 
limitations of many-body methods. As one approaches the particle emission thresholds, 
it becomes increasingly important to describe correctly the coupling to the continuum of 
decays and scattering channels [20l|2T]. While the full treatment of antisymmetrization 
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and short-range correlations has become routine in ab initio approaches to nuclear bound 
states, the many-body problem becomes more difficult when long-range correlations and 
continuum effects are considered. 

This paper is organized as follows. Section [2]outlines modern theoretical approaches 
to the nuclear many-body problem, starting with renormalization schemes for nuclear 
forces and ending with a discussion of various many-body methods, in particular aimed 
at describing systems near the nucleon drip lines. Section |3] presents selected results for 
stable and weakly bound systems, with an emphasis on long isotopic chains. Conclusions 
and suggestions for future work are contained in Sec. HI 

2. Theoretical foundations 

The task of developing a comprehensive theoretical framework that would be 
quantitative (i.e., capable of reproducing and explaining existing experimental data), 
have predictive power (i.e., capable of massive extrapolations), and provide uncertainty 
quantification (i.e., theoretical error bars) is daunting. To this end, theoretical models 
must meet three stringent requirements: (i) input (interactions, functionals) must be 
of high quality; (ii) many-body dynamics and correlations must be accounted for; and 
(iii) the associated formalism must take care of open-quantum-system aspects of the 
nucleus. 

Solving the nuclear many-body problem is a challenging endeavor, since "exact" 
solutions exists only for the very lightest systems and closed form solutions only 
for highly-idealized, non-realistic cases. Over the last few decades there has been 
tremendous progress in developing many-body ab initio techniques capable of treating 
light and medium-mass systems. The state-of-the-art methods are based on controlled 
approximations and the underlying computational schemes account for successive many- 
body corrections in a systematic way. 

The currently used theoretical techniques include: coupled-cluster methods [2214M] . 
Quantum Monte Carlo applications [251427] . perturbative expansions [28], Green's 
function methods [29l [30]. correlation operator methods [31], the density- matrix 
renormalization group [321431] . density functional theory [3 [35], in- medium similarity 
renormalization group [36|,[37], and large-scale diagonalization methods [38|[39]. 

In low-energy nuclear physics, baryons (protons, neutrons, and possibly deltas), and 
mesons (pions and possibly other mesons with masses below 1 GeV) are considered to be 
the appropriate degrees of freedom, meaning that one can derive effective Hamiltonians 
based on the interactions between these constituents. Those Hamiltonians are derived 
from Lagrangians that are consistent with the symmetries of low-energy QCD using 
chiral perturbation theory at different orders (z/), in terms of a perturbative expansion 
using a hard scale A and a soft scale Q [10]. Within this picture, three- nucleon forces 
(3NFs) arise naturally since the interactions are derived for effective degrees of freedom 
at low energy, and heavier baryons and mesons are integrated out. 

Currently, there are at least three different many-body methods for nuclear 
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structure with the capabihty to include both nucleon-nucleon interactions and 3NFs, 
yielding results that meet the few-body benchmarks. For light nuclei (with mass 
numbers A ~ 12), the Green's function Monte Carlo |40fH2] and large-basis no-core 
shell-model (NCSM) approaches [39lll3llll] have been successfully applied. These 
methods provide excellent results for two- and three-body Hamiltonians applied to 
light nuclei. However, for medium-mass and heavier nuclei, the dimensionality of the 
many-particle problem becomes intractable for these techniques. More recently, the 
coupled- cluster method has been applied to the structure of light and medium-heavy 
nuclei pl |H5] - H7] , providing excellent benchmarks for few-body systems such as ^He [23] . 
This technique is best suited for treating nuclei around closed-shells but has very 
advantageous scaling properties that enable accurate calculations in very large model 
spaces. At the singles and doubles level of the method (expanding Slater determinants in 
terms of the exponential of one-particle-one- hole and two-particle- two- hole correlations), 
the number of floating-point operation scales as n^n^, where hq and n„ are the numbers 
of hole and particle orbitals, respectively. Such soft scaling, when compared to the 
nearly combinatorial scaling of methods based on Hamiltonian diagonalization (as a 
function of basis size and/or particle number), allows one to build an extension of ab 
initio descriptions of nuclei all the way to medium and heavy systems. 

Nuclear ab initio methods are now evolving to tackle a crucial challenge; the 
description of open nuclear systems. This development is timely since the overarching 
scientific questions of modern nuclear structure are about very neutron-rich or proton- 
rich nuclei whose properties are impacted by a coupling to the particle continuum 
of scattering and decaying states. Many nuclei of interest lie every close to particle 
emission thresholds, i.e., they are either short-lived or unstable. Such nuclei cannot 
be described in a closed quantum formulation, which assumes that the nucleons are 
artificially confined in a trap; hence, they are unable to decay. 

It is therefore crucial for the nuclear many-body problem that new theoretical 
methods are developed that allow for an accurate description of loosely bound and 
unbound nuclear states. The recently developed complex-energy Gamow Shell Model 
(GSM) [38] has proven to be a reliable tool in the description of nuclei, where continuum 
effects cannot be neglected. In the GSM, a many-body basis is constructed from a 
single-particle Berggren ensemble [iQl - IST] . which treats bound, resonant, and scattering 
states on equal footing. Recently, GSM calculations of loosely bound and unbound 
states in nuclei, starting from a realistic interaction and a Gamow-Hartree-Fock basis 
were reported [52]. However, an ab initio description of light, unstable nuclei within 
the GSM approach will require novel many-body truncation schemes. In Ref. [53] the 
Berggren basis was employed for the first time in coupled-cluster calculations of the 
helium isotopes, and there are promising attempts to include the Berggren basis in an 
NCSM-like framework [53]. 

Finally, for the heavy and complex nuclei, the tool of choice is the nuclear Density 
Functional Theory (DFT) [55j and its extensions. Modern nuclear DFT is based on 
the mean-field approach rooted in the self-consistent Hartree-Fock-Bogoliubov (HFB) 
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problem. The DFT work is closely tied to ab initio studies of experimentally inaccessible 
systems such as neutron drops to enhance the predictive capabilities |56] - [58] . Important 
areas of research are the structure and decays of very neutron-rich nuclei, the dynamics 
of the fission process in heavy nuclei, and the structure of neutron star crust. 

The quasiparticle HFB energy spectrum contains discrete bound states, resonances, 
and nonresonant continuum states [59146T] . That is, the effect of the particle continuum 
is naturally incorporated into the formalism, provided that the system of interest is 
particle-bound. In order to treat weakly-bound nuclei accurately, special care should be 
paid to the spatial extension of HFB states [6l]. In this context, of particular interest are 
coordinate-space HFB approaches in large boxes [62], and PTG-HFB [63] and Gamow- 
HFB [6l] frameworks. 

The effective interaction in DFT, represented by the energy density functional 
(EDF), is characterized by a set of coupling constants constrained by experimental 
data and pseudo-observables determined by ab initio calculations. The uncertainty 
margins on EDF parameters are obtained by means of statistical methods like linear- 
regression with error analysis, which allows us to determine the correlations among EDF 
parameters, parameter uncertainties, and the errors of calculated observables [581I651468] . 
Such an approach is essential for providing predictive capability and extrapolability, 
and for estimating the theoretical uncertainties. A representative example of DFT 
calculations containing uncertainty quantification, highly relevant to this paper, are the 
large-scale DFT calculations of Ref. jS] assessing the limits of the nuclear landscape. 
Quantifying the limits of nuclear binding is important not only for understanding the 
mechanism of nuclear binding, but also for understanding the origin of elements in 
the universe. Indeed, the astrophysical processes responsible for the generation of many 
heavy elements operate very closely to the drip lines; hence, the structure of very exotic, 
weakly bound nuclei directly impacts the way the elements are produced in stars. 

In the following, we shall first briefly discuss models for the nuclear interactions and 
how these can be renormalized for truncated Hilbert spaces. Thereafter follows a brief 
description of configuration-interaction and coupled-cluster approaches to the nuclear 
many-body problem. We then discuss some developments that allow the study of open 
quantum systems within ab initio frameworks. 

2.1. Realistic nuclear interactions 

Ideally, the nuclear interactions should be derived from QCD, the underlying theory 
for the strong interaction. Unfortunately, the ^5*0 NN scattering length, a fundamental 
quantity characterizing the NN interaction, still remains a huge challenge for lattice 
QCD (LQCD) calculations, as the physical point (at the actual pion mass of 140 MeV) 
lies in a resonance in the unitary region [69]. In order to improve calculations, one 
needs to go to smaller lattice spacings, around b = 0.05 fm, but such computations 
require lattices of the order of 96'^ x 193 (the simulation cost grows as b~^). Meanwhile, 
the LQCD computations for NN systems [70] - [73] provide important insights about the 
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basic properties of the nuclear force. 

Several models of the NN interaction have been developed during recent years. 
These interactions reproduce NN scattering data up to 300-350 MeV laboratory energy 
with excellent precision [711476] . The number of free parameters used in the fitting 
procedure is around 40. Being optimized to on-shell scattering data, various interaction 
models introduce different off-shell behavior resulting in different predictions for finite 
nuclei. (For more discussion on off-shell effects, see, e.g., Ref. [77].) The interaction 
models of Refs. [711476] are either based on one-boson exchange with selected low-mass 
mesons, or are simply expressed in terms of local operators and optimized to data. More 
recently, due to a progress made in chiral-EFT [7H1I79] , several groups have constructed 
interaction models based on the underlying symmetries of QCD. 

The starting point of the chiral-EFT is an approximation to the QCD Lagrangian. 
The chiral effective Lagrangian is given by an infinite series of terms with an increasing 
number of derivatives and/or nucleon fields, with each term being dependent on the pion 
field according to the rules of the broken chiral symmetry. Applying this Lagrangian 
to the NN scattering results in a systematic series of Feynman diagrams in a small 
parameter, the ratio of Q/A^, where Q is pion mass and A-^^ ^ 1 GeV is the chiral 
symmetry breaking scale. For a given order z/, the number of contributing terms is 
finite and calculable; these terms are uniquely defined, and the prediction at each order 
is model independent. By going to higher orders, the amplitude can be calculated to 
any desired accuracy. This scheme has become known as chiral perturbation theory 
(XPT). 

Interactions derived from EFT have a number of advantages over traditional 
potential models. First, the corresponding currents are consistently formulated, and this 
is important for the correct description of observables other than the energy. Second, a 
power counting in terms of Q/A^ exists for systematic improvements of the interaction 
and observables. Third, the hierarchy of NN forces, 3NFs, and higher-rank forces is 
explained by power counting. Three-nucleon forces enter at next-to-next-to-leading 
order (N^LO) order, and four-nucleon forces appear at N^LO. Recent work [TOllSO] 
provides, for the first time, a chiral interaction of quantitative accuracy. The authors of 
Ref. [81] undertook the task of generating an accurate NN interaction based on chiral 
perturbation theory at N'^LO order. The number of free parameters used in this chiral 
interaction is 24, which is similar to the number of free parameters used to parametrize 
other NN forces. Three-nucleon forces have been explored recently in light nuclei and 
neutron and nuclear matter [82[[83]. In light nuclei, chiral 3NFs affect the binding 
energy, radii and transition probabilities [39]. They are responsible for the anomalous 
long half life of ^'^C [M]. In oxygen isotopes, they are believed to determine the position 
of the neutron drip line and the structure of neutron-rich isotopes [l5l[85]. In calcium 
isotopes, 3NFs are important for our understanding of shell evolution and (sub)shell 
closures [86],[87]. Three-nucleon forces also affect properties of neutron matter [881490] 
and the saturation of nuclear matter [83] . 



Living on the edge of stability, the limits of the nuclear landscape 



8 



2.2. Derivation of effective interactions for truncated Hilbert spaces 

The nucleon-nucleon interaction is strongly repulsive at short distances, meaning that a 
calculation which starts with the free interaction may converge very slowly. To overcome 
this problem, one defines effective interactions which can be used as starting points 
for calculations with basis-expansion approaches, such as the NCSM. The aim of this 
section is, therefore, to present and partly justify the computation of effective two-body 
Hamiltonians acting within a reduced Hilbert space. 

The starting point is a translationally invariant two- or three-body realistic 
interaction. The first method discussed is based on a diagonalization of the two- or 
three-body Schrodinger equation in a large harmonic oscillator (HO) basis. In the 
two-body case, this basis typically consists of several hundreds of HO shells. Via a 
similarity transformation, one obtains an effective two-body interaction, acting within a 
reduced model space used in NCSM. Alternatively, one can diagonahze the Schrodinger 
equation in a full momentum space and then carry out a similarity transformation to 
a smaller space defined by some momentum cutoff A. This leads to the so-called low- 
momentum approach, or l^ow-k, discussed below. We discuss also the so-called similarity 
renormalization group transformations. 



2.2.1. Similarity transformed effective interactions in an oscillator basis The 
translationally invariant Hamiltonian for an A-nucleon system, used in NCSM studies, 
is 



H 



E 



=1 



lPt_ 
2m 



2mA 



V^-Z^M is the three- nucleon force, and P 
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is the center-of-mass momentum. 
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Since we want to employ a HO single-particle basis, it is convenient to use the relation 

A 



mVt^ 
~2A 



0, 



(2) 



where R = l/A X]j=i is the center-of-mass coordinate and Vt is the oscillator frequency. 
Inserting this relation, we can rewrite the above Hamiltonian as H{Vt) = Hq + Hj — Hcm, 
where Hq is the HO Hamiltonian, Hj is the interaction term defined by 



mfl'^ 
~2A 



i<j<k 



ijk 

NNN^ 



(3) 



and Hqm = P'^/2mA + mAQ'^H'^ /2 is the center-of-mass term. 

Shell-model calculations are carried out in a model space defined by a projector P. 
The complementary space to the model space is defined by the projector Q = 1 — P. 
With the above Hamiltonian, we can then construct an effective interaction acting within 
the model space P, reproducing exactly Np eigenvalues of the full Hamiltonian. This 
can be accomplished by a similarity (Lee-Suzuki) transformation [91149^ . However, no 
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unique unitary transformation exists: one can construct infinitely many different unitary 
transformations wfiich decouple the model and complementary {Q = 1—P) subspaces, as 
discussed in Ref . |95] . Calculation of the exact A-hodj effective interaction is, however, 
as difficult as finding the full space solution. Using two-body interactions, the effective 
interaction is often approximated by a two-body effective interaction determined from 
a two-nucleon subsystem. 

The above-mentioned transformation can be performed also in a three-nucleon 
space. This will generate effective 3NFs even if the starting Hamiltonian includes only 
two-body terms. These effective interaction calculations are performed in a Jacobi 
coordinate HO basis |39l[96]. An effective three-body interaction can be computed 
also for starting Hamiltonians with pure NN interaction terms. This procedure has 
been shown to speed up convergence |96], although it should be noted that all cluster- 
approximated effective interactions (a < A) will in this approach reproduce the free 
interaction results in the infinite model space limit. 

2.2.2. Similarity transformed effective interactions in momentum space Alternatively, 
instead of carrying out the similarity transformation in an oscillator basis, one can 
perform the transformation in momentum space [97j. This approach consists of two 
steps: (i) a diagonalization of the two-body Schrodinger equation in the full momentum 
space, and (ii) a similarity transformation [93l|9l] to relative momenta k G [0,A] fm~^, 
with A defining the relative momenta model space. Typical values of A are in the range 
of ~ 2 fm~^. This evolution to lower /c- values can be performed using renormalization 
group (RG) methods [9H1[99]. We note, however, that the RG approach differs from the 
Lee-Suzuki transformation as the Q-space block of the effective Hamiltonian is now set to 
zero. The evolution to lower cutoffs shifts contributions from the sum over intermediate 
states to the interactions, just as RG equations in quantum field theory shift strength 
from loop integrals to coupling constants. The evolved low-momentum effective NN 
interactions have been dubbed Viow-k [97]. These interactions are significantly softer; 
hence, more perturbative, as seen, e.g., in nuclear matter calculations |100j . By 
construction, Viow-k preserves two-nucleon observables for relative momenta up to the 
cutoff. It should be noted, however, that no tractable approach to evolve many-body 
forces in this approach has been proposed so far. 

2.2.3. Similarity renormalization group transformed interactions In recent years, a 
new approach to perform the similarity transformation has been developed and used, 
in particular, with the chiral interactions. The similarity renormalization group (SRG) 
approach |101[|102j builds on the general principle of RG theory, namely that the relevant 
details of high-energy physics for calculating low- energy observables can be captured 
in scale-dependent coefficients of operators in a low-energy Hamiltonian. Using this 
principle as a guide, the transformed Hamiltonian Hs = UlHUs can be expressed 
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through a differential (ffow) equation 

^H...^J. (4) 

in which 7]s is the generator of the flow, an antihermitian operator that is related to the 
unitary transformation at the resolution scale s through 

- f - -4. (5) 

The evolved Hamiltonian is Hg = Tj-d + Vg, where the s-dependence is contained fully 
in the effective potential Vg. The crucial step of the method is the choice of the most 
appropriate generator 77. The antihermitian property is automatically fulfilled by 

Vs=[Gs,Hg], (6) 

where Gg is a momentum-diagonal operator. With this choice, the transformed operator 
flows towards a band-dagonal form in momentum space, thus achieving the desired 
decoupling property. Applications to nuclear forces have used Gg = Trci [97ill03 j . 

We should note that the generator of the flow is a many-body operator, which 
implies that it induces many-body terms in the effective operator, see Sec. 13. 1[ In 
principle, the evolution of many-body forces is straightforward, which is an advantage 
of the SRG method. However, the associated technical difficulties are significant. While 
the two-body evolution is usually performed in momentum space [103j . the three- 
body evolution was first implemented in a three-body HO basis |1041I105 ]. A first 
implementation of the SRG flow equation in three-body momentum space was first 
performed for a one-dimensional model [106j and, more recently, for realistic chiral 
interactions |107j . 

Finally, we want to highlight an important difference between the Lee-Suzuki 
transformation and the scale-dependent effective interactions discussed in Sees. 12.2.21 
and I2.2.3[ While the Lee-Suzuki transformation also induces many-body terms it will, 
by construction, reproduce the bare-interaction results already in a two-body cluster 
approximation (i.e., without induced many-body forces) - as the many-body space is 
extended towards the full Hilbert space. This property is not shared by the SRG and 
Viow-k effective interactions. By truncating the flow at two- or three-body level, the 
unitarity of the transformation is violated, and the bare result will not be reproduced, 
even when the many-body model space becomes the full Hilbert space. 



2.3. Configuration interaction and coupled cluster theory 

Present shell-model codes can reach dimensions of ~ 10^° basis states [38l[M fT08lfT09] . 
and Monte Carlo-based shell-model codes can attack problems with d ~ 10^^~^^ 
[26 |[271[TT0till2j . Although these numbers are impressive, the dimension limitation 
has important implications for calculations of nuclei that involve weakly bound states 
and/or resonances, as such states require still larger basis sets. Extensions of the NCSM 
to weakly bound systems will be discussed below. 
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Large-scale diagonalization is widely used in many areas of physics, from quantum 
chemistry |23] to nuclear physics |38]. The method is based on a projection of the model 
Hamiltonian onto a finite-dimensional subspace of the many-body Hilbert space; hence, 
the method is an instance of the Rayleigh-Ritz method. In the standard-shell model 
approach, one takes the stance that the many-body Hamiltonian is composed of two 
parts: Hq and Hj, treating the latter as a perturbation of the former. The eigenfunctions 
of the mean-field Hamiltonian, Hq, are assumed to comprise a single-particle basis for 
the Hilbert space. This leads to a matrix diagonalization problem, hence the name of 
the method. As Hj is the residual interaction that "perturbs" simple configurations of 
Ho, the method is commonly referred to as the configuration- interaction method. 

The NCSM method differs from the standard shell model in several ways: (i) there 
is no core, which implies that all nucleons are included in the many-body basis; (ii) no 
mean field is introduced and all interactions are included explicitly; (iii) the center-of- 
mass motion is separated exactly; and (iv) NCSM Hamiltonian contains realistic two- 
and three-body nuclear interactions and the Hamiltonian matrix is usually built from 
these interactions using the similarity transformations. 

Coupled cluster theory employs a different approach to systematically build the 
many-body wave functions using a large number of single-particle states. In actual 
calculations one typically limits the number of active many-body states to at most three- 
particle-three-hole excitations, whereas the single-particle basis can easily be extended 
to some 20 major oscillator shells. This has important consequences for studies of 
weakly bound systems. The coupled-cluster method has a rich history in both nuclear 
physics |113f[116j and quantum chemistry [23 tlll7j . 

In this contribution we highlight recent achievements and important lessons learned 
from ab initio calculations performed with the NCSM and coupled-cluster methods. In 
the following, we outline the essential features of these two approaches. 

2.3.1. No-core shell-model theory The starting point for NCSM calculations is the 
translationally invariant Hamiltonian H of Eq. ([T]). This Hamiltonian is however 
modified by adding a HO center-of-mass Hamiltonian Hqm = Tqm + Ucm, where Uqm = 
Amfl'^H'^ /2. In practice, the Hamiltonian H + Hqm is used as input to the computation 
of a similarity transformed effective interaction. At the stage of constructing the many- 
body Hamiltonian matrix, the HO center-of-mass term is subtracted and a Lawson 
projection term (3 {Hcm — f is added to shift spurious CM excitations up in the 
energy spectrum. 

The truncation of the many-body model space is usually performed in terms of the 
total energy. That is, the sum of HO excitations that is contained in a Slater determinant 
basis state is restricted according to X^f + ^i) < A^tot,max- More often, however, the 
basis truncation is specified in terms of iV^ax that measures the total number of HO 
excitations above the unperturbed ground state. For A = 3,4, these two measures are 
the same, but for p-shell systems they differ, e.g., for ^Li, A^max = Atot.max — 2, and for 
^^Li, A^max = Atot,max — 7, etc. As an example, the ^Li results that are presented in 
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Sec. 13.11 were obtained in model spaces up to A^max = 16. With this A^max-truncation, 
the highest single-particle state that can be reached has the energy {2n + l)hQ = 17hfl, 
which gives 171 possible nlj orbitals and 2280 individual single-particle states {nljm). 
The dimension of the many-body space is d = 0.79-10^. This Hamiltonian matrix is very 
sparse and can be diagonalized using iterative Lanczos algorithms for real, symmetric 
matrices |118] . The ground-state is usually reached with 10 — 20 Lanczos iterations. 
The method provides excitation spectra as well as wave functions, from which various 
observables can be computed. 

Since NCSM employs HO basis and the total energy truncation (including all 
allowed configurations), any eigenstate of the translationally invariant Hamiltonian 
factorizes into a product of a wave function depending on intrinsic coordinates, and 
a wave function depending on the center-of-mass coordinate. The use of any other 
single-particle basis, or any other truncation scheme, results in a mixing of intrinsic and 
center-of-mass motions. 

Several important ingredients are not considered here, such as the construction 
of antisymmetric, few-body states in the Jacobi-coordinate basis needed for the 
computation of interaction matrix elements. For the many-body systems considered 
here, the antisymmetric many-body states are constructed in the uncoupled M-scheme. 
For a much more complete description of the NCSM method, see Ref. [39] . 

2.3.2. Coupled cluster theory The single- reference coupled cluster theory is based on 
the exponential ansatz for the ground-state wave function of the A-nucleon system, 

l*o) = e^|$o), (7) 

where |$o) is an uncorrelated, closed shell reference state, and T is a linear expansion 
in particle-hole excitation amplitudes. 

The approximation that is made in the coupled-cluster approach is the truncation of 
the correlation operator T at a given (low-order) particle-hole excitation level. Note that 
in contrast to the full configuration interaction method, where the expansion in particle- 
hole excitation amplitudes is linear, the expansion is non-linear in the coupled-cluster 
approach due to the exponentiation of T. The most commonly used approximation 
in the coupled-cluster approach is the truncation of the operator T at the singles- 
and doubles excitation level (CCSD). Higher accuracy can be obtained by including 
triples excitations (CCSDT) [22 |I117] . In terms of computational cost, the CCSD 
method scales as n^n^, while the full CCSDT scales as nj^n^, where no represents the 
number of occupied orbitals and the number of unoccupied single-particle states. 
Since coupled cluster theory with inclusion of full triples CCSDT is usually considered 
to be too computationally expensive, several approximations to the solution of the 
CCSDT equations have been developed. A sophisticated way of approximating CCSDT 
is known as the A-CCSD(T) approach [11911120] . in which the left-eigenvector solution 
of the CCSD similarity-transformed Hamiltonian is utilized in the calculation of a non- 
iterative triples correction which scales as n^ra^. Excited states and neighbors of closed 
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shell nuclei can be computed within the equation-of-motion methods |117j . Recently, 
the coupled-cluster method in an angular momentum coupled scheme was developed 
in Ref. p ^ll21] . This allows to address medium mass and neutron rich nuclei starting 
from "bare" chiral interactions. Results presented in this contribution will focus on 
the CCSD and the A-CCSD(T) approach, with various equation-of-motion methods for 
open-shell nuclei and excited states p ^ll22] . 



2.4- Theoretical treatment of open systems 

As discussed earlier, the description of loosely bound and unbound nuclear states 
represents a challenge for nuclear structure models |123j . Several approaches to 
this problem are based on the use of the so-called Berggren basis [19145T] . Modern 
applications of the continuum shell model include the (real-energy) Shell Model 
Embedded in the Continuum |124H127] based on the Feshbach projection formalism 
[128] , and the (complex energy) GSM [l8lll29f[T3T| employing the Berggren basis 
The Berggren completeness relation can be written as: 



\Un){Un\ + I 



dk\uk){uk\ 



(8) 



where b are bound states, d - decaying resonant states, and the integral along a contour 
in the complex-fc plane represents the contribution from the non-resonant scattering 
continuum, see Fig. O In general, different contours can be used for different 
partial waves. 



Im(k) 



well-bound state 
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O bound 
O decaying 
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. ^ state 



(-\ narrow 
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Figure 2. Illustration of the Berggren completeness relation ([8]) in the complex k- 
plane. The bomid states are located on the positive imaginary axis. The weakly bound 
halo states lie close to the origin. The positive-energy resonant states are located in the 
fourth quadrant. Those with a small imaginary part can be interpreted as resonances. 
The complex-fc contour represents the non-resonant scattering continuum. 



The recent GSM applications include the analysis of threshold effects in 
multichannel coupling and spectroscopic factors |132[I133] . description of isospin mixing 
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in weakly bound nuclei |134] . explanation of behavior of charge radii in helium halo 
nuclei |135j . and studies of asymptotic normalization coefficients in mirror nuclei 
[136] . Another exciting development utilizing the Berggren basis is the extension of 
the coupled-cluster method to open systems; it represents the first complex-energy 
nuclear ab initio framework [l5],[53llHZ] capable of describing many-body bound and 
unbound states. In actual coupled-cluster calculations, one performs first a Hartree- 
Fock calculation and transforms the Hamiltonian to the Hartree-Fock basis. As a 
consequence, the coefficients t° of the singles amphtude, see Sec. 12.3.21 acquire small 
values in the solution of the coupled-cluster equation. One also finds, when using a 
HO basis, that this approach reduces the /if2-dependence of the computed energies. In 
practical computations one aims at increasing the number + 1 of employed oscillator 
shells until the results become virtually independent on the model-space parameters. 
In Ref. |52] it was demonstrated that one could start with a standard HO basis and 
corresponding interactions. We refer the reader to the latter reference for further details. 

In order to use the Berggren basis in large-scale configuration-mixing calculations, 
the integral over the non-resonant continuum is discretized by using a suitable 
quadrature rule. Using this discrete Berggren basis, the GSM basis is obtained in 
the usual way by constructing many-body Slater determinants. The dimension of the 
Hamiltonian matrix grows rapidly with the number of discretized continuum states and 
the number of nucleons; hence, advanced numerical methods that can handle large non- 
Hermitian matrices must be used. In the context of the GSM, it has been shown that 
the Density Matrix Renormalization Group (DMRG) is an efficient way to compute the 
low- lying spectrum of the Hamiltonian at a low computational cost. 

2.4-1 ■ Density Matrix Renormalization Group for the GSM The DMRG method was 
first introduced to overcome the limitations of the Wilson-type renormalization group 
to describe strongly correlated systems with short-range interactions [32 |ll37j . More 
recently, the DMRG has been reformulated and applied to finite Fermi systems [138] . 
nuclear shell model |139H141] . and open systems jl42j . While most of the DMRG 
studies have been focused on properties in strongly correlated closed quantum systems 
characterized by Hermitian density matrices, systems involving non-Hermitian and non- 
symmetric density matrices can also be treated [3 ^I142[[H3] . 

Let us consider the application of the J-scheme DMRG in the context of the 
GSM (GSM+DMRG). The objective is to calculate an eigenstate jj'^) of the GSM 
Hamiltonian H with angular momentum J and parity vr. As \J'^) is a many-body 
pole of the scattering matrix of H, the contribution from non-resonant scattering shells 
along the continuum contour to the many-body wave function is usually smaller 
than the contribution from the resonant orbits [H]. Based on this observation, the 
following separation is usually performed |142] : the many-body states constructed from 
the single-particle poles form a subspace A (the so-called 'reference subspace'), and the 
remaining states containing contributions from non-resonant shells form a complement 
subspace B (see Fig. [3]). 
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Figure 3. Schematic illustration of the GSM+DMRG procedure during the s step 
of the warm-up phase. States {fc^} from A, previously optimized states a^, and states 
{{lj)s} constructed by occupying the s*'' shell with n particles are coupled to generate 
the new set of states {fc^ ® is}''' = {kA ® {as ® (Ij)^}}'^- 

One begins by constructing states forming the reference subspace A. All 
possible matrix elements of suboperators of the GSM Hamiltonian H acting in A, 
expressed in the second quantization form, are then calculated and stored and the 
GSM Hamiltonian is diagonalized in the reference space to provide the zeroth-order 
approximation |\E'j)'^'') to \J'^). This vector, called 'reference state', plays an important 
role in the GSM+DMRG truncation algorithm. The scattering shells (Ij), belonging 
to the discretized contour L+, are then gradually added to the reference subspace to 
create the subspace B. This first stage of the GSM+DMRG procedure is referred 
to as the 'warm-up phase. For each new shell that is added, all possible many-body 
states denoted as \i)B are constructed and matrix elements of suboperators of the GSM 
Hamiltonian acting on \i)B are computed. By coupling states in A with the states 
lis), one constructs the set of states of a given . This ensemble serves as a basis in 
which the GSM Hamiltonian is diagonalized. The target state |\E'j) is selected among 
the eigenstates of H as the one having the largest overlap with the reference vector 
|\E'j)'^''). Then, the desired truncation is performed in B by introducing the reduced 
density matrix, constructed by summing over the reference subspace A |144] . In standard 
DMRG applications for Hermitian problems, where the eigenvalues of the density matrix 
are real non-negative, only the eigenvectors corresponding to the largest eigenvalues are 
kept during the DMRG process. Within the metric defining the Berggren ensemble, the 
GSM density matrix is complex-symmetric and its eigenvalues are, in general, complex. 
As a consequence, the truncation is done by keeping the eigenstates (the 'optimized' 
states) with the largest nonzero moduli of eigenvalue [M] . 

The warm-up phase is followed by the so-called sweeping phase, in which, starting 
from the last scattering shell (/j)iast, the procedure continues in the reverse direction 



Living on the edge of stability, the limits of the nuclear landscape 



16 



(the 'sweep-down' phase) until the first scattering shell is reached. The procedure is 
then reversed and a sweep in the upward direction (the 'sweep up' phase) begins. The 
sweeping sequences continue until convergence for target eigenvalue is achieved. 

Applications of GSM+DMRG have been reported in Refs. [311 [1351 [Tl2]. Those 
examples demonstrate that weakly-bound and unbound nuclear systems, which - 
because of a prohibitively large size of Fock space, cannot be treated by means of 
direct diagonalization techniques - can be treated very efficiently using DMRG. The 
DMRG also opens up the avenue to perform NCSM calculations for open systems using 
the Berggren basis |54j . 

2.4-2. The Resonating Group Method An alternative formulation for describing open 
quantum systems, characterized by a limited number of open channels, is the resonating 
group method (RGM) [145] . Here, the many-body wave function is decomposed into 
contributions from various channels that are distinguished by their different arrangement 
of the nucleons into clusters. In principle, this corresponds to the expansion of the 
bound state, or the interior region of a scattering state, into an over complete set of 
basis functions. The basis functions consist of two parts: the cluster wave functions and 
the wave function representing the relative motion of the clusters. In the case of two 
clusters, for instance, the full wave function of the system can be written as 

= 5^ A{$i.$2.¥^.(r-;)}, (9) 

V 

where v labels cluster channels. This expansion is complicated due to the presence 
of the antisymmetrizer Ay^ which accounts for the exchange of nucleons between the 
clusters. The intrinsic wave functions, $1,2, are internally antisymmetric and would, 
in the NCSM+RGM approach [HHIITIT], be eigenstates of the NCSM Hamiltonian for 
that particular cluster. 

At this point, a basis of binary-cluster channel states, |$;^"^), is introduced that 
includes the spin-coupled product of internal (antisymmetric) wave functions of the two 
clusters at relative distance r. This basis can be used to expand the many-body wave 
function: 

l^j^T^ = ^1 drr-'^-^^A, \^i7) . (10) 

By diagonalizing the NCSM Hamiltonian in this basis, one obtains the RGM equations: 
Y^ldrr^ {ni:j{r\T)-EMt:j{r\T)\ = 0, (11) 

where M'^/l,'^ {r' , r) and 'H'l7{r',r) are the norm and Hamiltonian kernels, respectively. 
These non-local quantities contain all the nuclear structure information. However, 
we note that the basis states are asymptotically orthogonal so that all important 
physical quantities (such as the scattering matrix) can be defined with the asymptotic 
solution. The non-orthogonality mainly appears at short distances and is primarily 



Living on the edge of stability, the limits of the nuclear landscape 



17 



due to antisymmetrization effects. The RGM equations can be solved by means of the 
standard symmetric orthogonahzation method |147] . 

The NCSM+RGM approach has been initially used for single-nucleon projectiles, 
both for scattering and bound-state cases |146m48j . The method has recently been 
extended to study projectile-target binary-cluster states, where the projectile is a 
deuteron |149j . This new development makes it possible to compute, e.g., d+'^B.e 
scattering, clusterized ^Li states, as well as deuteron fusion reactions on and ^He [150] . 

3. Selected results 

In this section, we focus on selected results for several chains of neutron-rich isotopes, 
both light and heavy. We base our theoretical analyses on large-scale NCSM 
calculations, coupled-cluster calculations, and nuclear DFT. An important question 
is the role played by three or more complicated many-body terms in the low-energy 
effective nuclear Hamiltonian. The study of 3NFs in medium mass and neutron rich 
nuclei is an ongoing research topic in nuclear many-body physics. Several calculations 
indicate that 3NFs are needed to understand various properties of nuclei, see for example 
Refs. [SHlllSlllHllinilHnilHlEHllHTlITnS]. As one moves to more neutron-rich isotopes, 
correlations play an increasingly more important role. The mean field contribution is 
often reduced relatively to the effects due to two-, three- and many-body correlations. 
Furthermore, close to the drip lines, the number of degrees of freedom increases 
dramatically, with resonant and non-resonant continuum channels becoming available. 
All of this requires a good understanding of correlations and interactions driving the 
observed properties. 

A three-body force is expected to play a role in the evolution of single-particle 
energies as more and more particles are added. For instance, the effect of the monopole 
term has been analyzed intensively in terms of phenomenological interactions over the 
last two decades [38lll51ffl53] . However, there is also clear evidence from several 
calculations that a traditional shell-model picture with a strong mean field that defines 
a single-particle basis, may break down, see for example Refs. [2T11H7] and discussion 
below. 

As outlined in this contribution, we are now in a position where 3NFs from EFT 
can be included routinely into various ab initio methods. This means that we can start 
to explore, which components of the nuclear forces are at play when we move towards 
the drip lines. Furthermore, with long chains of isotopes to be studied theoretically and 
experimentally, one can eventually attempt to extract in-medium information that will 
allow us to constrain 3NFs for heavier nuclei. 

In the following, selected results for light and medium-heavy are discussed. We 
shall focus on results for binding energies and a few excited states for selected chains 
of isotopes In several cases we will study the impact of 3NFs and couplings to the 
continuum. 
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3.1. Results for light nuclei 

In this section, we mainly focus on ab initio description of light nuclei near the 
drip lines. We shall restrict discussion to approaches that try to solve the many- 
body Schrodinger equation without uncontrolled approximations, using a well-defined 
microscopic Hamiltonian as a single input. While there are several realistic nuclear 
interactions that reproduce NN scattering phase shifts and few-body data with high 
precision, we focus mainly on results obtained with chiral interactions. 

Bound state ab initio calculations of light nuclei were pioneered by the Variational 
Monte Carlo (VMC) and Green's Function Monte Carlo (GFMC) methods [25]. They 
have studied A = A — 12 systems using various combinations of Argonne NN interactions 
and UIX or Illinois (IL) 3NFs, and have demonstrated an impressive agreement with 
experimental energies, not only for bound states but also for narrow resonances |154j . 

Ab initio methods that can currently model systems with A > 4, include 
the Effective Interaction for Hyperspherical Harmonics (EIHH) [1551 1156] . lattice 
calculations [8011157] . the coupled-cluster method [MllllGllTSSj . and the No-Core Shell 
Model [391 ill [159]. The EIHH is currently limited to A < 6 due to the difficulty of 
antisymmetrization. Calculations on the lattice for light bound systems [157] and the 
Hoyle state [80] have been recently reported. 

As a benchmark example, and an illustration of the use of modern chiral interactions 
including consistent 3NFs in NCSM, we present in Table [1] results for A = 3,4 obtained 

Table 1. Properties of "^H and ^He. Benchmarking of calculations with chiral 
NN+3NF interactions using NCSM |160] and Hyperspherical Harmonics jl61] methods. 



(From Refs. [39l[l62].) 





NN (N^LO) 


+3NF(N2L0) 


Expt. 




NCSM HH 


NCSM HH 




3H [MeV] 


7.852(5) 7.854 


8.473(5) 8.474 


8.482 


(r2)V2 [fm] 


1.650(5) 1.655 


1.608(5) 1.611 


1.60 


^He Eg, [MeV] 


25.39(1) 25.38 


28.34(2) 28.36 


28.296 


(r2)V2 [fm] 


1.515(2) 1.518 


1.475(2) 1.476 


1.467(13) 



with and without inclusion of 3NFs. The chiral 3NF is from order N^LO and contains 
two low-energy constants (LEC) that need to be determined from A > 2 data. In 
Ref. [160] . these constants were obtained from A = 3 binding energies, but also by 
investigating sensitivities of properties of A> 4 nuclei to the variation of the constrained 
LECs. Besides the triton ground state energy, which is by construction within a few 
keV of experiment, the NN+3NF results for the '^He ground-state energy and point- 
proton radius are in excellent agreement with measurement. We also note a very good 
agreement between NCSM and the variational HH method [161] . 

Developments based on an importance-truncation (IT) have recently been proposed 
in Ref. [163] to address the factorial growth of the NCSM model space with A/'max and 
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particle number A. With a criterion based on perturbation theory, a large fraction of 
the many-body basis states can be discarded and calculations with A^max = 12 for ^^C 
and ^^O have become feasible [105] . 

Let us first consider the chain of He isotopes. The charge radii and masses of He 
isotopes (up to ^He) were recently determined experimentally [16411165] and compared 
to results from ab initio methods. For ^He in particular, results for radii and binding 
energies are available from calculations with GFMC [I66], NCSM [I67], FMD [T68] . 
and EIHH [16511169] . Unfortunately, at the present stage, those results can not be used 
for benchmarking as different interactions have been used. Actually, the GFMC results 
are the only published converged calculations that include 3NFs. Within the range of 
GFMC results, it is possible to reproduce both the charge radius and separation energy. 
However, there is a large uncertainty due to different models of 3NFs and different trial 
wave functions used. 

The first ab initio coupled-cluster calculations using a Berggren basis were 
performed in Ref. [170] . Within the coupled cluster singles-and-doubles approximation 
(CCSD), the ground state binding energies and lifetimes of the ^"^"He isotopes were 
calculated employing Viow-k with cutoff A = 1.9fm~^, derived from the N^LO nucleon- 
nucleon interaction of Ref. [8T]. A comparison between the CCSD results and experiment 
is shown in Fig. [Hand demonstrates a fair agreement. An improved description of the 
helium isotopes must include the effects of 3NFs and triples correlations (CCSDT). 
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Figure 4. Left panel: CCSD ground state energies (red solid line) compared to 
experimental binding energies (black dashed line) for '^"^'^He. Right panel: Imaginary 
parts of the CCSD (red solid line) and experimental (black dashed line) ground state 
energies for'^^^^He. 



A hallmark feature of halo systems appearing near the nucleon drip line is the rapid 
increase of matter radii. Corresponding effects on the charge radius are much smaller, 
but this quantity can be studied very accurately using laser spectroscopy techniques. 
Let us consider the neutron-rich nuclei ^'^He as an example. The charge radius of ®He 
is 2.059(7) fm and exceeds that of ®He, which is 1.959(16) fm, see Refs. [TMfT65] . Both 
of them are much larger than the charge radius of ^He, 1.681(4) fm. These two heavy 
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helium isotopes are neutron halo nuclei and, since the protons are confined inside the 
tightly bound a-core, the differences of charge radii carry unique structural information 
on the nuclear Hamiltonian and many-body dynamics. In Ref. |135] . ^He and ^He were 
described with the GSM as systems of valence neutrons moving outside the a core, 
and interacting via a finite-range Minnesota potential. The core-neutron potential was 
described by a Woods-Saxon potential and the model space was constructed from the 
Berggren basis. While such a calculation is obviously not not ab initio, it allows us 
to understand intricate experimental data in simple terms. Specifically, the GSM work 
demonstrated that the observed charge radii depend mainly on three factors: (i) the 
recoil due to the motion of valence nucleons around the a-core; (ii) the spin-orbit term; 
and (iii) the swelling of the a-core in the neutron environment (see Fig. [5]). While (i) 
and (ii) are robustly predicted by GSM (i.e., they depend weakly on the interaction, 
provided that the threshold energies are under control), the core swelling effect (iii) 
must be taken from an ab initio theory (here, from GFMC calculations [166] ). One can 
thus conclude that the GSM approach links the high-quality atomic data to a subtle 
in-medium effect that provides a stringent test of ab initio theory. 




Figure 5. Different contributions to the charge radius of ^He (solid line, squares) and 
^He (dashed line, dots) calculated in GSM. The core swelling contribution is taken 
from GFMC calculations of Ref. [53]. Experimental charge radii come from |165] 
(triangles). The NCSM 167_ (stars) and GFMC |166) (pentagons) results are marked 
for comparison. The inset shows GSM rms neutron radii compared to experiment [171) . 
(From Ref. [HS].) 



The chain of lithium isotopes offers many splendid examples of drip line physics. 
The nucleus ^Li belongs to the valley of beta-stability and ^'^'^Li are particle-bound but 
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beta-unstable. Adding one more neutron to ^Li yields an unbound ^°Li, while adding a 
pair of neutrons produces the additional pairing energy that makes ^^Li bound. Such 
an even-odd staggering effect in binding energy can be often seen along the neutron 
drip line. The conventional shell model picture of ^^Li puts the valence neutrons in a 
p-wave state. However, measured angular correlations in fragmentation reactions [172] 
demonstrate the presence of states with different parities. Furthermore, the structure 
of the unbound ^°Li is expected to correspond to a virtual s-state with a very large 
scattering length. These phenomena: melting and re-organizing of shell structure, 
ground states embedded in the continuum, and dilute matter densities such as halos, 
appear in most isotopes around the drip lines. They just illustrate a point that will 
be made clear throughout this contribution, namely, light nuclei living on the edge of 
stability cannot be described within a mean-field picture. 

The Li and Be isotopic chains were recently studied in the NCSM \173\ to investigate 
the wealth of exotic properties that generally pose a challenge for nuclear-structure 
models: (i) appearance of clustering; (ii) halo structure of ^^Li; and (iii) low ground-state 
quadrupole moment of ^Li. In this systematic study, series of calculations were carried 
out in large model spaces using CD-Bonn 2000 |74j , and INOY (IS-M) |174j interactions. 
The degree of convergence was estimated from the A^max- and dependence of the 
results. Alternatively, one can utilize the fact that NCSM calculations performed at 
different HO frequencies should all converge to the same value in the limit Amax — ^ oo. 
This constitutes an example of multiple converging sequences in the NCSM, discussed 
extensively in Ref. jl75j . In this context, we note that much work is currently focused 
on quantifying the basis truncation corrections [1 ^1176] . 

Figure [6] shows the ground-state energies of A = 6 — 11 Li isotopes. The isotopic 
trend is nicely reproduced, but also the known feature of pure NN interactions giving too 
little binding for many-nucleon systems is seen. The INOY interaction is a bit different 
as the NN P-wave scattering has been modified slightly in order to reproduce binding 
energies and the analyzing powers in the A = 3 systems. Note that the model spaces 
used in computations of ^^Li were not large enough to reach the exponential convergence 
region. 

In Fig. [7] we compare the calculated and experimental trends for a number of 
observables for the Li chain of isotopes. With the important exception of the radius of 
the ^^Li halo ground-state, a good agreement between NCSM results and experiment 
is found. In particular, the overall trends are well reproduced. In this context, 
however, it is important to point out that calculations of electromagnetic moments are 
usually performed in the impulse approximation (i.e., one-body nucleon currents), and 
that the inclusion of meson exchange currents can add about 20% to magnetic dipole 
moments [177] . We expect more work in this direction, including fully consistent two- 
body currents from chiral Lagrangians. Another success of NCSM calculations, is the 
reproduction of the small quadrupole moment of ®Li that is known to pose a challenge 
for theory. In particular, the general failure of three-body models for this observable has 
been blamed on missing antisymmetrization of the valence nucleons and the nucleons in 
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Figure 6. Ground-state energies for Li isotopes predicted in NCSM compared with 
experiment. The exponential convergence rate was not fully reached for -'^^Li. (From 
Refs. [TO^ITTg] .) 




Figure 7. Electric quadrupole moments, magnetic dipole moments, and charge radii of 
Li isotopes obtained in NCSM and compared with experiment. (From Refs. [16211173) .) 



the alpha-core |178j . In addition, the ratio Q (^^Li) /Q (^Li) is found to be very close to 
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unity, as confirmed by very precise experimental data |179] , and the trend for tlie mucli 
larger moments of A = 7 — 11 is nicely reproduced. This has been obtained within the 
truncated HO basis space that does not give a very accurate description of the dilute 
halo structure of ^^Li. Still, the observed decrease of the charge radius in A = 6 — 9 is 
reproduced. 

A clear example of the disappearance of magic numbers in weakly bound systems is 
found in ^^Be. The experimental ground state of ^^Be is an intruder 1/2+ level, while the 
first p-shell state is l/2~ situated at = 320 keV. The neutron separation energy is only 
503 keV, and there are no other bound states. An investigation of A = 9, 11 isotopes in 
large-basis ab initio NCSM calculations were reported in Ref. |180] . Calculations were 
performed for both natural-parity and unnatural-parity states in model spaces up to 
A^max = 9 using four different accurate NN potentials. The ^^B Amax = 9 calculation, 
with a basis dimension of 1.1 ■ 10^, was the largest NCSM diagonalization at that time. 
The calculations did not reproduce the anomalous 1/2+ ground state, but did predict a 
dramatic drop in the positive-parity excitation energies with an increasing model space. 
Furthermore, the behavior of the INOY results suggested that a realistic NNN force 
might have an influence on the observed parity inversion. 

It is expected that the halo character of a loosely bound neutron should have a 
major influence on the characteristics of a nuclear state. This fact is illustrated by 
the extremely strong El transition between the two bound states in ^^Be, that can be 
explained by a large overlap of the initial and final state wave functions at large distances. 
The El strength was underestimated by a factor of 20 in the NCSM calculations |180] , 
demonstrating that the halo character of this state is extremely hard to reproduce 
using a HO basis. This observation makes ^^Be an excellent candidate for testing the 
NCSM-I-RGM method, in which the relative motion of the core and the valence nucleon 
is treated more accurately. By imposing bound-state boundary conditions to the set 
of coupled channel Schrodinger equations of Eq. ( ITTj) . the bound states of ^^Be could 
be studied in a NCSM-I-RGM model space spanned by the n+^^Be channel states with 
inclusion of the ground state plus three excited states of ^°Be |146] . Binding energies of 
the various ^°'"Be states are displayed in Table [2] for both NCSM and NCSM-^RGM. 
The precise treatment of the neutron halo strongly influences the 5'-wave relative kinetic 
and potential energies. The rescaling of the relative wave function in the internal region, 
seen in the NCSM-I-RGM, is the main cause of the dramatic decrease (~ 3.5 MeV) of 
the energy of the 1/2+ state. This effect makes the 1/2+ state bound and even leads to 
a parity inversion in NCSM-I-RGM. 

Let us mention also modern nuclear reaction calculations that offer opportunities to 
compute new observables that probe other properties of realistic nuclear Hamiltonians. 
Here we present results from A^ -|- a scattering using the NCSM-I-RGM approach. The 
A = 5 system is an ideal testing ground for many-body scattering theory for several 
reasons: (i) the A = 5 system does not have a bound state; (ii) "^He is tightly bound 
so that single-channel scattering is valid up to ~ 20 MeV; (iii) there are two low-lying 
p-wave resonances (3/2^ and 1/2^); (iv) non-resonant s-wave scattering (1/2+) for 



on the edge of stability, 


the 


limits of the 


nuclear 


landscape 












"Be(l/2-) 




1^ ) 


-^max 


E 


E 


Eth 


E 


Eth 


NCSM |180| 


8/9 


-57.06 


-56.95 


0.11 


-54.26 


2.80 


NCSM [180J 


6/7 


-57.17 


-57.51 


-0.34 


-54.39 


2.78 


NCSM+RGM |146| 






-57.59 


-0.42 


-57.85 


-0.68 


Expt. 




-64.98 


-65.16 


-0.18 


-65.48 


-0.50 



Table 2. Energies (in MeV) of the ground state of ^"^Be and the lowest states of 
"Be, calculated in NCSM using the CD-Bonn NN potential [74] at Srj = 13 MeV. The 
NCSM+RGM results were obtained using n+^'^Be configurations with A^max = 6 g.s., 
2+ 2+, and l| states of i°Be. (From Refs. [1461118(1] . 'l 
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Figure 8. Calculated phase shifts for (left panel) n-a and (right panel) p-a scattering, 
using the N'^LO NN potential [HI], compared to an i?-matrix analysis of data (+). 
Theoretical results include the ^He g.s., 0+0, 0^0, 1^0, 2~0, and 2~1 states. 

(From Refs. [H^fT^ .l 



which large effects of the Pauh exclusion principle is expected; (v) it is a well studied 
system. In particular, there have been recent microscopic studies of -|- a scattering 
with GFMC [m] and NCSM+RGM [TMfW] . 

A comparison with an accurate i?-matrix analysis of the nucleon-a scattering is 
presented in Fig. |H1 It reveals that for both neutron (left panel) and proton (right panel) 
projectiles one can describe well the and, qualitatively, also the ^.03/2 phase shifts, 
usmg the N^LO NN potential. The good agreement of the N^LO ^5*1/2 phase shifts 
with the i?-matrix analysis can be credited to the repulsive action of the Pauli exclusion 
principle at short nucleon-a distances, which masks the short-range details of the nuclear 
interaction. On the other hand, the same interaction is not able to reproduce well the 
two P-wave phase shifts, which are both too small and too close to each other. 

Recently, the NCSM-(-RGM approach was applied to a low-energy radiative capture 
reaction '^Be(p, 7)®B |149] . and two fusion reactions ^H((i, n)^He and ^He((i,p)^He |150] . 
In all SRG evolved chiral EFT potential at two-body level was considered. It 

is important to realize that reaction calculations are extremely sensitive to the exact 
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positions of thresholds and resonances. In the studies mentioned above, the authors 
could utihze the model-dependence on the two-body SRG flow parameter to reproduce 
the experimental separation energies. This will, however, remain a particular challenge 
for ab initio reaction calculations starting from nuclear Hamiltonians which do not 
involve free parameters. 

3.2. Neutron rich oxygen and calcium isotopes 

The oxygen isotopic chain has been extensively studied experimentally during the last 
years, see for example Refs. |182[ll83j . The neutron drip line has been tentatively 
established at ^^O, but the particle stabihty of ^^O is still a matter of a debate since the 
current experiments are on the limit of the production cross section. The isotopes ^^O 
(A^ = 14) and ^^O (A^ = 16) exhibit a doubly magic nature. Their structure is believed 
to be governed by the evolution of the lsi/2 and 0^5/2 orbits. Recent experiments show 
that 25,26q g^j^g unbound with respect to ^^O |1841I185] . On the other hand, the two- 
neutron drip line for the fluorine isotopes extends beyond ^^F. One speculates that the 
behavior of the neutron drip line for the oxygen and fluorine isotopic chains arise from 
a delicate balance between the proton-neutron and neutron-neutron interactions, the 
coupling to the continuum and 3NFs [45 | [85 | l87]. 

Another chain of isotopes crucial for theoretical developments is the calcium chain. 
It contains several possible closed-shell nuclei beyond the well established ones, namely 
^°Ca and ^^Ca. The A^ = 32 sub-shell closure has been established from experiments 
on calcium |186[I187] . titanium |188] . and chromium |189j . The nucleus ^^Ca has a 
reduced value of 2f excitation (but more than twice as large as seen in open-shell 
calcium isotopes) than that observed in ^^Ca, suggesting a sub-shell closure. For 
^^Ca (A^ = 34), no sub-shell closure has been seen experimentally in chromium |190] 
or titanium |191lll92j . and there are some doubts regarding a sub-shell closure in 
calcium |193] . The heaviest calcium isotopes that have been observed are ^^'^^Ca |194] . 
but the masses have been measured only up to ^^Ca |195] . 

As emphasized in Sec. [21 to describe physics of nuclei near the drip lines, one 
needs to properly take into account (i) the effects of 3NFs, (ii) the presence of open 
decay channels and particle continuum, and (iii) many-nucleon correlations. Recently 
coupled-cluster calculations have been carried out for the binding energies and spectra 
of the neutron rich oxygen and calcium isotopes taking for the first time all these effects 
into account. The effects of 3NFs were included effectively in terms of density dependent 
corrections to the NN interaction. Those were derived from 3NFs by summing over the 
third particle in symmetric nuclear matter. This is obviously a departure from a rigorous 
ab initio approach, but nevertheless it is a first step towards including effects from 
both 3NFs and coupling to the particle continuum in coupled-cluster calculations. The 
continuum effects were included by using a Berggren basis for the relevant partial waves. 
The binding energies per particle for selected oxygen and calcium isotopes are displayed 
in Fig. |9l The inclusion of 3NFs significantly improves agreement with experiment. In 




Figure 9. Binding energy per nucleon for selected oxygen and calcium isotopes 
calculated in coupled cluster theory with NN (triangles) and NN+SNF (red squares) 
forces, compared to experimental data (circles). (Based on Refs. [151 157] .) 



order to investigate the evolution of shell structure and binding energy systematics in 
the neutron rich calcium isotopes in more detail, we show in Fig. [10] the total binding 
energies for very neutron rich calcium isotopes ranging from ^'^Ca to ^^Ca. Several 
interesting features are apparent: (i) 3NFs are important; (ii) a peninsula of weak 
binding appears for the very neutron rich calcium isotopes ^°'^^'^^Ca, and (iii) the results 
with NN interactions only, do not yield a flattening of binding energy at large neutron 
numbers. Looking in more detail at the structure of ^°^^^Ca, we find the J'^ = 1/2^ 
ground state of ^^Ca slightly above the threshold, unbound by about 0.2 MeV with 
respect to ^''Ca, and entirely dominated by an s wave. Furthermore, the ordering of the 
gds shells is found to be reversed as compared to the standard shell model ordering. In 
particular, we find that the 5/2^ state appears below the 9/2'*' state in ^^'^^'^^Ca. It is 
to be noted that the calculations employing the oscillator basis yield for ^^Ca the level 
ordering of the conventional shell model, with a 9/2'*' ground state spin assignment, 
thus suggesting strong effects due to the continuum coupling. The nuclei ^^Ca and ^^Ca 
are predicted to be only weakly unbound with respect to ^°Ca. 

Figure [11] shows the results for the 2'^ states in the neutron rich calcium isotopes. 
It is interesting to note that coupled-cluster calculations predict only a weak sub-shell 
closure in ^^Ca with a 2"*" state around 1.9 MeV. These examples clearly point tp the 
need for 3NFs in neutron-rich nuclei. Overall, we find that the effects of 3NFs and 
the scattering continuum are essential for understanding the evolution of shell structure 
towards the drip line. Although more investigations are clearly needed, our results hint 
at a situation where odd calcium isotopes beyond ^°Ca are unbound, while even isotopes 
can be weakly bound. 

It is interesting to relate the coupled-cluster results for the Ca isotopes to those 
based on nuclear DFT. Figure [12] shows the two-neutron separation energies for even- 
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Figure 10. Binding energies of the neutron-rich calcium isotopes calculated in 
coupled cluster theory with NN (diamonds) and NN+SNF (squares) forces, compared 
to experimental data (circles). The dashed line indicates the NN result normalizes at 
48 Ca. (Based on Ref. [87].) 




Figure 11. First 2+ excited state energies for the calcium isotopes 42,48,50,52,54,56(^g^_ 
Black circles: experimental data; blue diamonds: results from nucleon-nucleon 
interactions; red squares: results including the effects of three-nucleon forces. (Based 
on Ref. [57].) 
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even calcium isotopes computed in Ref. [S] using SLy4, SV-min, UNEDFO, and UNEDFl 
EDFs and obtained in FRDM [I96] and HFB-21 |197] mass models. All those models 
predict consistently the neutron drip line around ™Ca. Interestingly, in all cases the 
two-neutron drip lins (5'2n = 0) is approached fairly gradually in all cases; this indicates 
that the neutron chemical potential A„ ~ — S'2n/2 stays close to, but below, zero for 
60 < < 70. A similar result was obtained in other DFT calculations jl9811201] . As 
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Figure 12. Theoretical extrapolations towards drip lines for the two- neutron 
separation energies S2n for the isotopic chain of even-even calcium isotopes using 
different EDFs (SLy4, SV-min, UNEDFO, UNEDFl), and FRDM [196] and HFB- 
21 |197) mass models (see text for more details). All models are consistent with the 
available experimental data. Detailed predictions around S2n = are illustrated in 
the inset. Drawn by Erik Olsen based on Ref. [S]. 



discussed in Refs. [6T |l202fl205j . the persistant appearance of A„ just below the threshold 
can be associated with the continuum effect due to pairing. Indeed, the scattering of 
neutron pairs into the close-lying non-resonant continuum gives rise to a stabilization of 
binding energy and a very weak dependence of A„ on A^; hence, extension of the range of 
bound nuclei. In terms of HFB quasiparticles, this continuum coupling manifests itself 
through increased occupations of low-lying quasi-particle states, especially those having 
low orbital angular momentum |202j . as one approaches the threshold. This results in 
a gradual increase of contribution from nonresonant continuum to the ground state of 
the system. 

By looking at canonical HFB states, one can notice the emergence of bound 
canonical orbits from the single-neutron continuum. For the considered case of the 
drip line Ca nuclei, the canonical neutron states Si/2, ^5/2, and 6/3/2 appear very close to 
the bound gg/2 level, which is expected to form the valence shell in the traditional shell 
model |198] . The fact that high-£ and low-£ orbits bunch up very close to the A„ = 
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threshold creates an opportunity for other correlations to further lower the binding 
energy. Figure [13] shows the isoscalar (mass) ground-state quadrupole deformations 
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Figure 13. Isoscalar quadrupole deformation /32 for Z < 50 even-even nuclei predicted 
in HFB calculations with UNEDFO energy density functional [67]. The neutron-rich 
Ca isotopes with iV sa 50 are predicted to be deformed. (Based on Ref. [8].) 

in Z < 50 even-even nuclei calculated in the HFB-I-UNEDFO model in Ref. [8]. The 
appearance of deformation around ™Ca. is clearly seen. For more discussion on deformed 
drip line nuclei, see Refs. p [205ti208] . 

The binding energy stabilization close to the threshold appears naturally in the 
continuum shell model through the collective coupling of shell model states via the 
decay channel |209[I210] . This coupling, governed by the anti-Hermitian term in the 
effective Hamiltonian that represents the continuum coupling, leads to the formation of 
the collective aligned state. The mechanism responsible for the creation of an aligned 
state is similar to the formation mechanism of super-radiant states |211[I212] . In this 
language, the behavior of drip line Ca isotopes, the particle stability of ^'^He and 
^^Li, and the appearance of cluster states near the reaction threshold that exhausts 
most of the decay width, are all manifestations of the same emergent near-threshold 
phenomenon |209[I210| . 

Clearly, it would be very interesting to see if ab initio approaches confirm the trends 
predicted in nuclear DFT. The fact that the coupled-cluster calculations of Ref. [87] yield 
a J"^ = 1/2+ ground state for ^^Ca, unlike the shell ordering in the conventional shell 
model, is a tantalizing hint that this is indeed the case. 
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4. Conclusions and perspectives 

In this contribution, we have demonstrated that modern many-body techniques are 
nowadays capable of providing a reliable description of nuclear properties thanks to new 
conceptual insights as well as algorithmic and computational advances. For light nuclei, 
like the chain of lithium and beryllium isotopes, no-core shell-model calculations provide 
an invaluable tool to link nuclear structure to the underlying forces. Similarly, coupled 
cluster theory that includes 3NFs and coupling to the particle continuum, is capable of 
providing reliable predictions for heavier nuclei such as oxygen and calcium isotopes. 
Modern nuclear density functional theory, with optimized energy density functionals, 
provides an excellent description of heavy nuclei, and it links to ab initio approaches 
as well. Recent parametrizations of nuclear energy density functionals provide fairly 
consistent predictions when extrapolated to mass regions where experimental data are 
not available. Close to the limits of stability, the degrees of freedom represented by 
resonances, weakly bound states and the non-resonant continuum need to be accounted 
for properly. The latter has important consequences for the interpretation of rare 
isotopes in terms of a naive shell-model picture. Our results indicate that this traditional 
picture may not be relevant close to the neutron drip line. All of this holds great promise 
for a quantitative and predictive modeling of nuclei from a bottom-up perspective. 

It is important to emphasize that the nuclear many-body problem is a splendid 
example of a multiscale problem, with length scales spanning many decades. A 
description of multiscale processes entails different theoretical methods for different 
length scales. First-principle methods are limited to few interacting nucleons. With 
an increasing number of degrees of freedom, DFT-based methods become the methods 
of choice. For even larger systems, a molecular-dynamics-based modeling is the favored 
approach. To link these different scales and methods properly is a great challenge not 
only to nuclear physics, but to fields as diverse as material science and life science, see 
for example Ref. |213j . 

In this work, we have illustrated some challenges and opportunities of the modern 
nuclear many-body problem, especially in the context of rare isotopes. In particular, 
we have emphasized the need for a high quality input, the importance of many-body 
dynamics and the impact of the coupling to open channels, summarized in Fig. [HI 
With a fundamental picture of nuclei based on the correct microphysics, we can remove 
the empiricism inherent today, giving us thereby a greater confidence in the science we 
deliver and the predictions we make. Guided by unique data on rare isotopes, we are 
embarking on a comprehensive study of all nuclei, based on the most accurate knowledge 
of nuclear interactions, the most reliable theoretical approaches, and massive use of new 
computer hardware and advanced numerical algorithms. The prospects are excellent. 
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Figure 14. Challenges of the nuclear many-body problem. A comprehensive 
theoretical framework that would be quantitative, have predictive power, and provide 
uncertainty quantification must meet three stringent requirements: (i) the input must 
be quantified and of high quality; (ii) many-body dynamics and correlations must 
be accounted for; and (iii) the associated formalism must take care of open-quantum- 
system aspects of the nucleus. Only then can we hope to understand rare isotopes, such 
as ^^Li (two-neutron halo), ^^O (doubly-magic, probably unbound), ^''Ca (territory for 
new physics, where ab initio theory and DFT meet), and ^ssy (r-process system, fission 
recycling participant). 
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